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The Self-Consistent Ornstein-Zernike Approximation (SCOZA) is an accurate Uq- 
uid state theory. So far it has been tied to interactions composed of hard core 
repulsion and long-range attraction, whereas real molecules have soft core repulsion 
at short distances. In the present work, this is taken into account through the in- 
troduction of an effective hard core with a diameter that depends upon temperature 
only. It is found that the contribution to the configurational internal energy due to 
the repulsive reference fluid is of prime importance and must be included in the ther- 
modynamic self-consistency requirement on which SCOZA is based. An approximate 
but accurate evaluation of this contribution relies on the virial theorem to gauge the 
amplitude of the pair distribution function close to the molecular surface. Finally, 
the SCOZA equation is transformed by which the problem is reformulated in terms of 
the usual SCOZA with fixed hard core reference system and temperature-dependent 
interaction. 



I. INTRODUCTION 



The Self- Consistent Ornstein-Zernike Approximation (SCOZA) has been found to give 
very accurate results for the equation of state for fluids and lattice gases. This approximation 
was proposed by H0ye and Stell^'^. Some preliminary results were obtained for supercritical 
temperatures^. 

However, for subcritical temperatures there were numerical problems connected to the 
singular behavior along the spinodal curve and the no-solution region inside it. The first 
successful subcritical solution was made by Dickman and Stell for the Ising model or the 
lattice gase case^. It turned out that results were very accurate. 

Borge and H0ye made a more general numerical investigation of the SCOZA equation of 
state in the critical region^. Clearly scaling was not fulfilled very close to the critical point. 
But apart from that the critical behavior was close to that of real fluids. Especially it was 
noted that the critical exponent (3 for the curve of coexistence was equal or close to 0.35. 
This value was subsequently confirmed by numerical evaluations by others^. 

Also Borge and H0ye compared SCOZA results with experimental data on CO2 in the 
critical region'^. By closer study of these data on a scaling plot it was realized that they did 
not collapse onto a single scaling curve but gave a series of close-lying curves instead, and 
the SCOZA results were consistent with such deviations from scaling^. 

H0ye, Pini and Stell then made a closer investigation of the critical region of SCOZA^. 
Their analysis showed that SCOZA fulfills a generalized kind of scaling instead of the usual 
one, and they obtained the value 0.35 for the critical exponent {3. 

SCOZA was then applied to continuum fluids, and Pini, Stell and H0ye considered the 
hard sphere fluid with interaction of Yukawa form^. The results aggreed well with simulation 
data. Furthermore, interactions using more Yukawa terms were also considered^^. This 
allowed approximations to the Lennard- Jones (LJ) potential of realistic fluids to be studiedi^. 
Accurate results were also obtained for the more short-ranged fullerene interaction o^^'^^ . 

So far SCOZA has been applied to continuum systems with hard cores perturbed by an 
attractive interaction only. Real molecules, on the other hand, have soft cores. In order 
to obtain more accurate results, this should therefore be taken into account. To do so we 
use hard spheres with an effective diameter d. In combination with the self-consistency 
requirement of SCOZA, however, the use of an effective diameter turns out to become a 



non-trivial problem. This is connected to the necessity of properly taking into account the 
contribution to the internal energy due to the repulsive interaction. Furthermore, a nu- 
merical procedure using the semi-analytical solution of the Ornstein-Zernike (OZ) equation 
becomes more problematic as a grid of fixed densities p implies varying packing fractions 
7] = ^pd^ when the effective diameter d changes. While a fully numerical solution of the 
OZ relation is possiblo^^ and sidesteps this particular problem, it entails a huge increase in 
computational complexity and cost. 

In the present work, on the other hand, we avoid the need for heavy numerics. Instead we 
opt for a more conceptual approach to the introduction of soft cores into SCOZA by means 
of effective hard cores. In doing so we prefer here a prescription for the hard core diameter d 
that depends only upon temperature ( sect ion HT]). It turns out that an accurate evaluation of 
the reference system contribution to the configurational internal energy plays an important 
role. This quantity can be linked to the pressure and hence to the virial theorem f section UTT]) 
and must be taken into account for the energy route part of the SCOZA self-consistency 
problem (section llV|l . A further refinement uses the virial theorem once more, this time 
to gauge the amphtude of the pair distribution function close to the core (sectional)- We 
conclude with a formal mapping of the soft-core SCOZA for fixed interaction onto the 
usual SCOZA with fixed hard core and a temperature dependent potential. This essentially 
eliminates any softness and so provides a way of avoiding the complexities associated with 
it (section EH). 



II. INTERACTION AND EFFECTIVE HARD CORES 



The pair interaction for real atoms that are neutral is typically given by the LJ interaction 
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with length and energy scales of a and e, respectively. (From now on we will use units where 
cr = e = 1.) In approximations, the repulsive part of this interaction where 0(r) > is 
commonly replaced by a hard core. This has also been done in earlier SCOZA computations 
where the hard core diameter d was kept fixed. 

State-dependent effective hard cores have been used earlier to evaluate the equilibrium 
properties of fluids, and there exist various recipes; see Ref. 3 for a recent compilation and 



comparison in the context of thermodynamic perturbation theory. The most simple schemes 
such as the Barker- Henderson prescriptioniS use an effective diameter d that depends on 
temperature only, whereas more advanced recipes include a density dependence, too^. 

In the present work we limit ourselves to a c? that depends only upon temperature. The 
main reason for this is the added complexity of the problem when combining it with the self- 
consistency of SCOZA. In particular, the latter requires that the reference system internal 
energy becomes part of the thermodynamic self-consistency and must be properly accounted 
for. This inclusion of the contribution of the soft cores in the self-consistency problem is 
the fundamental difference between this method and another one which is a combination of 
SCOZA and a first order perturbation theory proposed as a way of handling soft repulsive 
cores^I and studied recently^^. 

Another reason for the use of a c? that depends only upon temperature is the connection 
to the hard sphere equation of state that becomes unique. In other words, both the virial 
theorem and the compressibility theorem will remain consistent if they are consistent for 
the reference system with temperature independent d. 

The importance of d for quantitative accuracy is seen from the fact that the mean-field 
critical temperature and density vary with d as d^ and d~^ respectively when the interaction 
ip{r) is held constant. From the second virial coefficient of the repulsive part of the LJ 
interaction one will find roughly d^ ~ 0.92 at the critical temperature if 0(r) for r < 1 
is replaced by an effective hard core. This is in satisfactory agreement with the SCOZA 
result with fixed d when compared with simulations. With temperature independent d = 1 
SCOZA yields the value Tg = 1.245 with a 3- Yukawa fit to the Lennard- Jones interaction for 
r > 1, cf. Fig. 4.3 of Ref. |l9l For simulations performed for the full LJ interaction, one has 
found Tc = 1.31020 and 1.3132^ by molecular dynamics, and = 1.3120(7)22 and 1.326(2)^ 
by Monte Carlo methods. Previous SCOZA evaluations but with temperature independent 
d = 1, however, yielded values Tc = 1.304, 1.293^^1, and 1.305^ of the critical temperature. 
But these values are based on 2- Yukawa fits to the LJ interaction with a compensation for 
the soft core such that the right second virial coefficient is reproduced. The variations in 
these latter numbers is due to the precise potential fit and the detailed form of the direct 
correlation function used with respect to the reference system hard cores. 

Without thermodynamic self-consistency, results for isotherms at different temperatures 
are independent. The effective d for a given temperature then determines the equation of 



state for that temperature. This is the case for common fluid theories. For SCOZA, on the 
other hand, this is no longer true. But away from the critical region this coupling cannot 
be important as other theories are accurate there, too. One approximate way to implement 
SCOZA with temperature-dependent d can then be to solve the equations with the same d at 
all temperatures but to use the results only for the isotherm corresponding to the chosen d. 
Repeated evaluations with different d will then give the full phase diagram. This approach 
bears some resemblance to the SCOZA-based perturbation theory mentioned above^^. Just 
as the latter, however, such a procedure is only approximate, and compiling a full phase 
diagram requires substantial computer resources far in excess of those needed for a single 
run of the SCOZA program. 

The division of the pair interaction 0(r) into a reference system part (t)o{r) and a per- 
turbing part il){r) can be performed in various ways. As suggested above, one can use the 
repulsive part where 0(r) > as the reference system interactiooi^. 




00 (r) = 

With this the perturbing interaction is 

i,{r) = { - (1) 




However, there are more choices for this splitting of the interaction, and an alternative is to 
use2i (rm = 2^) 

0(r) - 0(rm) r < 
r > 1 



ilj{r) 



()(r) r > Tin. 

As was verified by H0ye and Borge, SCOZA requires a perturbing interaction that is mainly 
attractive, or else the equations cannot be solved^^. Clearly, the above suggested splittings 
of the potential fulfill this condition. 

Furthermore, in order to keep the numerical implementation less demanding it is desirable 
to approximate the interaction outside the effective hard core as a sum of Yukawa terms. 
This allows the OZ equation to be solved in a semi-analytic way. For a potential like the 



LJ one, such a multi- Yukawa form of 0(r) outside the core is easily found. For example, a 
simple non-linear least-squares fit of a sum of three Yukawa terms constrained to reproduce 
0(1) = converges rapidly and gives a result that is essentially indistinguishable from the 
original LJ form for r > 1^^. When d < 1, however, this fit is certainly not constant for 
d < r < 1 as mandated by the two types of splitting mentioned above. For instance one 
can then add one more Yukawa term to approximate the desired form; this will generally be 
a rather short-ranged function that hardly contributes beyond r = 1. For prescription (^J 
both the range and the amplitude of the additional Yukawa term can be fixed by imposing, 
e. g., ip{d) = ip'{d) = 0. At any rate, any remainder of ilj{r) not accounted for by this sum 
of Yukawa terms is added to (poi^) = 4'{r) — "ipir) and so enters the computation through 
the evaluation of the effective diameter. 

In solving SCOZA, the reference system is used as a boundary condition at temperature 
T oo, or /5 = where /3 = l/^ksT) and ks is Boltzmann's constant. Strictly speaking, 
the reference system becomes the ideal gas in this limit for a soft core, i. e., d ^ which is 
far from unity. However, use of the ideal gas as reference system may give rise to additional 
problems, especially in the numerical implementation. So we have not tried to investigate 
this possibility further. Instead we have focused upon the situation with an effective hard 
core diameter d near 1. The justification for this lies in our arguments above. They mean 
that different temperatures and densities do not couple significantly away from the critical 
point anyway. We can therefore start at /3 = with a non-vanishing value of d corresponding 
to, say, twice the critical temperature, and start to vary d only at lower temperatures for 
which SCOZA values then will be valid. 



III. REFERENCE SYSTEM AND REPULSIVE INTERNAL ENERGY 

In approximating a soft repulsive core by an effective hard core, the simplest prescription 
is to let the diameter d depend only upon temperature. One possibility is to define d such 
that the second virial coefficient for the soft particles coincides with that of the effective 
hard cores. 




The prescription of Ref. is slightly different from this in that it replaces the three- 
dimensional integral Q with a one- dimensional one. However, the result will be the same 



to leading order in the difference 1 — d which is considered small. In this connection it can 
be mentioned that the precise prescription is not crucial since both give a d that depends 
only upon /5, and the SCOZA problem needs only d{l3) as input, not its prescription. 

At high density there are better approximations with density dependence, but we expect 
the deviations of such choices from Eq. ^ to be small when the soft potential can be 
considered steep at the molecular surface. 

In contrast to strict hard spheres, a soft repulsive potential implies that the reference 
system also contributes to the configurational internal energy. This can be related to the 
pressure and ultimately to the virial theorem. The equation of state can be written as 

Pp = P+(i-p4-] lo, (3) 
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where p is the pressure, Iq = —(3pfo, and /o is the excess (beyond the ideal gas) Helmholtz 
free energy per particle. For effective hard cores Iq/p only depends on the packing fraction 
V, 

h = py{v)- 

Insertion into Eq. Q thus gives 

(3p = p(l-^^). (4) 



drj ^ 

The configurational or excess internal energy uq per particle of the reference system is now 

dio dy drj 

dy d In d^ 
^ ~^d^ 8(3 

and finally, using Eq. (jU for dy/dr], 

d In d? 

puo = (f^P-p)—gp-- 
The virial theorem for hard spheres implies 

271 

PP- P = —d^p'^nQ{d+) 

and so connects the above internal energy to the contact value nQ{d+) of the pair distribution 
function, 

pUo = yfl P no{d+) . (5) 



To gain some intuition for the consistency of this result, we can insert the definition 
for the effective diameter. This yields 

dlnd^ _ 1 dd^ _ J 0o(r)e-^*«Wdf 
dp ~ d^W ~ 

and thus finally 

puo = lp'no{d+) j 0o(r)e-^<^°(^)dr. 

This is the low density value of the internal energy, amplified by the contact value nQ{d+) 
for higher densities. Compared with the exact Mq, the pair distribution function is here 
approximated by rio((i+)e~'^"^"*-'"^ The accuracy of this increases as the repulsive part of the 
potential becomes less soft. A different prescription for d[(3) such as that of Ref. Q gives 
different expressions for the above two integrals. But otherwise the precise choice of d{l3) is 
of no consequence for the remainder of this work. 



IV. SCOZA EQUATIONS 

The SCOZA approach is based upon thermodynamic consistency between the energy and 
compressibility routes to thermodynamics. These routes are connected via the thermody- 
namic relation 



da d'^put 
where 



P^^' (6) 
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is the reduced inverse compressibility and Ut is the total configurational internal energy per 
particle. Both a and Ut are evaluated from the pair correlation function in different ways, viz., 
by the compressibility and energy routes to thermodynamics. In general, the pair structure 
is known only approximately, and a and Ut give different thermodynamics. In SCOZA, on 
the other hand, consistency between the two routes in the form of Eq. (jH)) is enforced by 
adjusting an unknown parameter, usually the amplitude of the direct correlation function 
c(r) outside the hard core. Specifically, in the Mean Spherical Approximation (MSA) the 
contribution to c(r) from the perturbing attractive interaction ipijr) is 



c(r) = —(3ip{r), r > d. 



(7) 



SCOZA replaces (3 in the above relation by an effective value that depends on both tem- 
perature and density and is obtained from the solution of Eq. The total correlation 
function h{r) = n{r) — 1 is then obtained from the core condition 



h{r) 



-1, 



r < d, 



and the OZ equation 



h{k) =d{k) + p5{k)h{k). 



(The tilde marks Fourier transforms.) One can then evaluate a as well as u, the internal 
energy contribution from the attractive interaction alone, according to 



In this way, both a and u are in principle obtained as functions of the effective temperature 
/3e for any given density and interaction. Consequently, u (or a) can replace /3e as the free 
parameter so that correspondingly a (or u) can be expressed in terms of the former. Inserting 
this into Eq. (jH)) one then obtains the SCOZA partial differential equation (PDE) for u (or 



V. CONTACT VALUE OF THE PAIR DISTRIBUTION FUNCTION 

For soft repulsive interactions the total internal energy ut to be used in Eq. (jnj is the 
sum of the contributions (0) and (jH)), i- e.. 



At first sight one might expect Uq to play a minor role. But a closer investigation shows 
that it is crucial for obtaining results consistent with a changing diameter d. The reason for 
this is that the reference system enters primarily as the boundary condition of the PDE at 
(3 = 0. For (3 > 0, a change in d is not "seen" except through uq (and a small perturbation 
of u due to small changes in h{r)). Neglecting uq thus means essentially keeping d fixed at 
its /5 = value. Our preliminary numerical work strongly indicated this importance of the 
Uq term. 



a = 1 — pc(0). 




(8) 



a). 



Ut = Uo + u. 



(9) 



Since Mq plays such an important role, the accuracy of the contact value n{d+) = l+h{d+) 
becomes of interest. For j3 > 0, n{d+) deviates from nQ{d+), and clearly the former more 
accurately describes the energy due to the reference system interaction (polr). With a direct 
correlation function of the form (0), however, there is no reason to expect that the h{r) 
obtained is accurate close to the hard core. In particular, the contact value h{d+) strongly 
depends on the choice of iplr) a.t r = d and thus on the potential fit with the added Yukawa 
tails. A more reliable method of obtaining the contact value is thus desirable. 

One appealing possibility is provided by the virial theorem that then also, to a certain 
degree, enters the SCOZA where it has played no role traditionally. H0ye and Stell earlier 
proposed full consistency between the energy, virial and compressibility routes-, but this 
requires structure functions depending on two free parameters and has not been considered 
numerically so far, nor will it be considered here. But on the other hand we can still obtain 
desired information about the contact value n{d+) via the virial theorem using the SCOZA 
quantities as input. 

From the virial theorem for hard cores with attraction, the contact value n{d+) follows 

as 

Ait 

PP = P+ —d^p^n{d+) + pp{u - v), 
o 

where the virial integral is split into two parts, with u given by Eq. (jH}. The expression for 
the remainder v then becomes 

If 

V = u -\ — p / (fVip(r))n(r)dr 
6 J 

1 



p / V (fijj(r))n(r)df, 
6 J 

with the integrations restricted to r > d. Note that v vanishes in the mean field limit as 
then n{r) = 1. Expression (jSj) for uq can now be expressed in terms of the new contact value 
n{d+) instead of the reference system value no{d+). With the above relations we find 

d In d"^ 

puo = {(3p- p- (3pu + (3pv) . (10) 
Thus the additional complication of evaluating uq is the evaluation of the integral for v. 



VI. TRANSFORMATION TO UNIT DIAMETER 



With varying d there is an additional problem if the analytic solution of the OZ equation 
for a sum of Yukawa terms is used. The latter provides a as a function of u only at constant 



packing fraction and interaction whereas the PDE requires the temperature derivative to 
be taken at constant density. When (3 and thus d change while the density grid in the 
discretization is kept fixed, the packing fraction also changes and the relation between a and 
u becomes less direct and more challenging to evaluate numerically. 

If now the varying d problem can be transformed into a situation of fixed unit diameter 
d = 1, the discrepancy mentioned will not arise and application of the analytic solution of 
the OZ equation will be simpler. The price to pay is an interaction that varies with j3 in the 
transformed problem. As we will see below, this works out nicely and gives equations that 
can be given a direct physical interpretation. 

To obtain the desired transformation we introduce a number of quantities, marking those 
of the unit diameter problem by a subscript 1: 

pi = pd^, Pi = p/d^, pi = pd\ 
Ui = ud^, Vi = vd^, (t'lij'i) = d^'t'ij')^ 
Ti = r/d, dri = dr/d^, ni(ri) = n(r). 

With Eqs. ©, ®, and (HHI) we have the SCOZA PDE 



da d'^ 
dp^^dp 



where for brevity we put 

dlnd^ _ p dd^ 

~^~d(r ~ d^w 

For a d that depends only on P, the introduction of the new quantities on the right hand 
side is done in a straightforward way. With d"^ /dp^ = d^d'^/dpf we get 



da 1 d 
dp " d^^^Jpl 



For a we find likewise 

dm diPiPi] 



dp dpi 
The derivative with respect to P then becomes 

da dPi dai ^ dpi dai 



ai. 



dp dp dpi dp dpi 

( \ P dd^ \ dai j3 1 4 ^'^1 

d^'d^'dp) Wi^^ P^ dfh 

,.dai ^ d'^pi 
I - A)— + piA- 



1 

d^ 



dPi dpi 



Upon insertion into Eq. (fTT|) the pi terms cancel and we obtain 

(9a 1 (9^ 



Pi — [Piui + Aipivi] 



where 



The latter equality follows from 
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dd^ 



d%l + A,). 



df3i " ' ' d(3i 

The resulting Eq. (|T2|) can be given a direct physical interpretation in terms of the 
transformed system. The latter consists of hard spheres of fixed diameter {di = 1), where 
now there is an attractive pair interaction that depends upon the temperature, = 
d^(f){r) = d^(f){rid). Furthermore, the soft repulsive interaction is no longer present as the 
transformed system has hard cores. In the usual virial graph expansion, the Helmholtz free 
energy per particle at given density will depend only upon the interaction as before (besides 
temperature and density). This will remain the same even though 0i(ri) depends on the 
temperature. The total configurational internal energy per particle obtained from the pair 
distribution function is thus 

5(/?i/i) 

d 
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d/3i 



Pi 
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[(3id^(j){rid)] ni(ri)dfi 

1 
3 



0i(ri) +Ai(0i 



0i(ri) + -AiV(fi0i(ri)) 



dd 
ni{ri)dfi 



ni{ri)dfi 



= Ui + AiVi, 
in accordance with the right hand side of Eq p2p . 



VII. CONCLUSION 



We have investigated a method to perform SCOZA evaluations for realistic molecules with 
soft cores. A temperature dependent effective hard core diameter is then introduced. For 
thermodynamic self-consistency it turns out that the excess internal energy of the reference 
system is important and should be treated accurately. For this purpose the contact value 



of the pair correlation function at the hard sphere surface is needed. Since SCOZA does 
not give such a rehable contact value, it is instead obtained by use of the virial theorem. 
In section IVII the SCOZA problem with varying d is transformed to a simpler one with 
fixed d = 1. This transformation not only eliminates some of the difficulties associated with 
temperature-dependent diameter but also provides an independent method of evaluation 
that can be useful as a test of the internal consistency of results. 

The results of this work have been obtained along with numerical work to implement and 
solve new problems when soft cores are considered. So far we have only considered simple 
test functions d{(3) to obtain a program that can handle varying d when using analytic 
expressions. In this respect the transformation to fixed = 1 in section IVTl has been verified 
numerically by solving SCOZA with both the original and transformed equations by putting 
for simplicity the virial type integral for v and thus the one for vi equal to zero. For explicit 
evaluations with the LJ interaction it is necessary to evaluate integral © for d{l3) and 
the less trivial one for v\ preliminary computations already including v are encouraging. 
Furthermore, it will then also be desirable to study the infiuence of different prescriptions 
for d{(3) and for the splitting of the potential into a reference part and a perturbing attractive 
part. We intend to extend our computations in this way. 
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